A modified ziggurat algorithm for generating exponentially- and normally-distributed pseudorandom numbers

The Ziggurat Algorithm is a very fast rejection sampling method for generating PseudoRandom Numbers (PRNs) from statistical distributions. In the algorithm, rectangular sampling domains are layered on top of each other (resembling a ziggurat) to encapsulate the desired probability density function. Random values within these layers are sampled and then returned if they lie beneath the graph of the probability density function. Here, we present an implementation where ziggurat layers reside completely beneath the probability density function, thereby eliminating the need for any rejection test within the ziggurat layers. In the new algorithm, small overhanging segments of probability density remain to the right of each ziggurat layer, which can be efficiently sampled with triangularly-shaped sampling domains. Median runtimes of the new algorithm for exponential and normal variates is reduced to 58% and 53% respectively (collective range: 41-93%). An accessible C library, along with extensions into Python and MATLAB/Octave are provided.